home *** CD-ROM | disk | FTP | other *** search
/ Developer CD Series 2000 November: Tool Chest / Dev.CD Nov 00 TC Disk 1.toast / Sample Code / Devices and Hardware / Velocity Engine / VelEng FFT / Test vBigDSP.c < prev    next >
Encoding:
C/C++ Source or Header  |  2000-09-28  |  17.9 KB  |  726 lines  |  [TEXT/CWIE]
  1. /*
  2.     Disclaimer:    IMPORTANT:  This Apple software is supplied to you by Apple Computer, Inc.
  3.                 ("Apple") in consideration of your agreement to the following terms, and your
  4.                 use, installation, modification or redistribution of this Apple software
  5.                 constitutes acceptance of these terms.  If you do not agree with these terms,
  6.                 please do not use, install, modify or redistribute this Apple software.
  7.  
  8.                 In consideration of your agreement to abide by the following terms, and subject
  9.                 to these terms, Apple grants you a personal, non-exclusive license, under Apple’s
  10.                 copyrights in this original Apple software (the "Apple Software"), to use,
  11.                 reproduce, modify and redistribute the Apple Software, with or without
  12.                 modifications, in source and/or binary forms; provided that if you redistribute
  13.                 the Apple Software in its entirety and without modifications, you must retain
  14.                 this notice and the following text and disclaimers in all such redistributions of
  15.                 the Apple Software.  Neither the name, trademarks, service marks or logos of
  16.                 Apple Computer, Inc. may be used to endorse or promote products derived from the
  17.                 Apple Software without specific prior written permission from Apple.  Except as
  18.                 expressly stated in this notice, no other rights or licenses, express or implied,
  19.                 are granted by Apple herein, including but not limited to any patent rights that
  20.                 may be infringed by your derivative works or by other works in which the Apple
  21.                 Software may be incorporated.
  22.  
  23.                 The Apple Software is provided by Apple on an "AS IS" basis.  APPLE MAKES NO
  24.                 WARRANTIES, EXPRESS OR IMPLIED, INCLUDING WITHOUT LIMITATION THE IMPLIED
  25.                 WARRANTIES OF NON-INFRINGEMENT, MERCHANTABILITY AND FITNESS FOR A PARTICULAR
  26.                 PURPOSE, REGARDING THE APPLE SOFTWARE OR ITS USE AND OPERATION ALONE OR IN
  27.                 COMBINATION WITH YOUR PRODUCTS.
  28.  
  29.                 IN NO EVENT SHALL APPLE BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL OR
  30.                 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
  31.                 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  32.                 ARISING IN ANY WAY OUT OF THE USE, REPRODUCTION, MODIFICATION AND/OR DISTRIBUTION
  33.                 OF THE APPLE SOFTWARE, HOWEVER CAUSED AND WHETHER UNDER THEORY OF CONTRACT, TORT
  34.                 (INCLUDING NEGLIGENCE), STRICT LIABILITY OR OTHERWISE, EVEN IF APPLE HAS BEEN
  35.                 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
  36. */
  37.  
  38. #include <MacMemory.h>
  39. #include <math.h>
  40. #include <stdio.h>
  41.  
  42. #include "vBigDSP.h"
  43.  
  44. #ifndef PI
  45. #define PI (3.1415926535897932384626433832795)
  46. #endif
  47.  
  48.  
  49. #pragma mark -
  50.  
  51. static double
  52. RMSE_real(
  53.     float *    x,
  54.     float *    y,
  55.     int         n
  56. )
  57. /* Calculates the root-mean-squared error between two real signals
  58.  * x[n] and y[n].
  59.  */
  60. {
  61.     double         t, err=0.0;
  62.     int         i;
  63.     
  64.     for (i=0; i<n; i++)
  65.     {
  66.         t = x[i] - y[i];
  67.         err += (t * t);
  68.     }
  69.     
  70.     return(sqrt(err/n));
  71. }
  72.  
  73.  
  74. static double
  75. RMSE_complex(
  76.     float *    x,
  77.     float *    y,
  78.     int     n
  79. )
  80. /* Calculates the root-mean-squared error between two complex signals
  81.  * x[n] and y[n].
  82.  */
  83. {
  84.     double         t, err=0.0;
  85.     int         i;
  86.     
  87.     for (i=0; i<n; i++)
  88.     {
  89.         t = x[2*i] - y[2*i];
  90.         err += t * t;
  91.         t = x[2*i+1] - y[2*i+1];
  92.         err += t * t;
  93.     }
  94.     
  95.     return(sqrt(err/n));
  96. }
  97.  
  98. #pragma mark -
  99.  
  100. static void
  101. ConvolveRealLiteral(
  102.     float *    x,
  103.     float *    y,
  104.     float *    z,
  105.     int         n
  106. )
  107. /* Performs a literal cyclic convolution on x and y, placing the
  108.  * result in z. This is an O(n^2) algorithm useful
  109.  * for accuracy testing.
  110.  */
  111. {
  112.     int            s, i, q;
  113.  
  114.     for (s = 0; s < n; s++)
  115.     {
  116.         z[s] = 0;
  117.         for (q = 0; q < n; q++)
  118.         {
  119.             i = (s-q)%n;
  120.             if(i<0) i+= n;
  121.             z[s] += y[i] * x[q];
  122.         }
  123.     }
  124. }
  125.  
  126. static void ConvolveComplexLiteral(
  127.     float *    x,
  128.     float *    y,
  129.     float *    z,
  130.     int         n
  131. )
  132. /* Performs a literal cyclic convolution on x and y, placing the
  133.  * result in z. This is an O(n^2) algorithm useful
  134.  * for accuracy testing.
  135.  */
  136. {
  137.     int            s, i, q;
  138.  
  139.     for (s = 0; s < n; s++)
  140.     {
  141.         z[2*s] = z[2*s+1] = 0;
  142.         for (q = 0; q < n; q++)
  143.         {
  144.             i = (s-q)%n;
  145.             if(i<0) i+= n;
  146.             z[2*s] += y[2*i] * x[2*q] - y[2*i+1] * x[2*q+1];
  147.             z[2*s+1] += y[2*i+1] * x[2*q] + y[2*i] * x[2*q+1];
  148.         }
  149.     }
  150. }
  151.  
  152. static void
  153. Convolve2DRealLiteral(
  154.     float *    x,
  155.     float *    y,
  156.     float *    z,
  157.     int         w,
  158.     int         h
  159. )
  160. /* Literal cyclic convolution in two dimensions, for testing purposes. The
  161.  * convolution of x[] and y[] is output in z[].
  162.  */
  163. {
  164.     int         k, j, q, p, kk, jj;
  165.     
  166.     for (k=0; k<h; k++)
  167.     {
  168.         for (j=0; j<w; j++)
  169.         {
  170.             z[j + k*w] = 0;
  171.             for (q = 0; q < h; q++)
  172.             {
  173.                 kk = (k-q<0)?k-q+h:k-q;
  174.                 for (p = 0; p < w; p++)
  175.                 {
  176.                     jj = (j-p<0)?j-p+w:j-p;
  177.                     z[j + k*w] += x[p + q*w] * y[jj + kk*w];
  178.                 }
  179.             }
  180.         }
  181.     }
  182. }
  183.  
  184. #pragma mark -
  185.  
  186. #define COMPLEX_TEST_MIN_POWER    0
  187. #define COMPLEX_TEST_MAX_POWER    18
  188.  
  189. static void TestComplexFFT()
  190. {
  191.     float     *real_data;
  192.     float     *vector_data;
  193.     int        i, j;
  194.     long    currentLength;
  195.     OSErr    result;
  196.     double    rmse_vector;
  197.  
  198.     for (i = COMPLEX_TEST_MIN_POWER; i<=COMPLEX_TEST_MAX_POWER; i++) {
  199.         
  200.         currentLength = 1 << i;
  201.  
  202.         real_data = (float*)NewPtr(2*currentLength*sizeof(float));
  203.         vector_data = (float*)NewPtr(2*currentLength*sizeof(float));
  204.         
  205.         if (real_data == nil || vector_data == nil) {
  206.             printf("\nerror allocating space for complex test");
  207.         } else {
  208.         
  209.             /////////////////////////////////////////////////
  210.             // initialize signal
  211.             /////////////////////////////////////////////////
  212.  
  213.             for (j=0; j<currentLength; j++) {
  214.                 real_data[2*j] = vector_data[2*j] =
  215.                     cos(PI*j*j/currentLength);
  216.                 real_data[2*j+1] = vector_data[2*j+1] =
  217.                     sin(PI*j*j/currentLength);
  218.                     
  219.             }
  220.             
  221.             /////////////////////////////////////////////////
  222.             // perform forward and inverse FFT and compare
  223.             // result with original signal
  224.             /////////////////////////////////////////////////
  225.             
  226.             result = FFTComplex(vector_data, currentLength, -1);
  227.             if (result != noErr) {
  228.                 printf("\nerror in forward complex FFT");
  229.             }
  230.  
  231.             result = FFTComplex(vector_data, currentLength, 1);
  232.             if (result != noErr) {
  233.                 printf("\nerror in inverse complex FFT");
  234.             }
  235.  
  236.             rmse_vector = RMSE_complex(real_data, vector_data, currentLength);
  237.             
  238.             printf("\n complex length %d FFT rmse: %g", currentLength, rmse_vector);
  239.             
  240.         }
  241.         
  242.         if (real_data) DisposePtr((Ptr)real_data);
  243.         if (vector_data) DisposePtr((Ptr)vector_data);
  244.     
  245.     }
  246.     
  247.     
  248.     
  249.     
  250. }
  251.  
  252. #define REAL_TEST_MIN_POWER    0
  253. #define REAL_TEST_MAX_POWER    18
  254. static void TestRealFFT()
  255. {
  256.     float     *real_data;
  257.     float     *vector_data;
  258.     int        i, j;
  259.     double    negativeOneToI = -1;
  260.     double    rmse_vector;
  261.     long    currentLength;
  262.     OSErr    result;
  263.     
  264.     for (i=REAL_TEST_MIN_POWER; i<=    REAL_TEST_MAX_POWER; i++) {
  265.         currentLength = 1 << i;
  266.         
  267.         real_data = (float*)NewPtr(currentLength*sizeof(float));
  268.         vector_data = (float*)NewPtr(currentLength*sizeof(float));
  269.         
  270.         if (real_data == nil || vector_data == nil) {
  271.  
  272.             printf("\nerror allocating space for complex test");
  273.  
  274.         } else {
  275.  
  276.             /////////////////////////////////////////////////
  277.             // initialize signal
  278.             /////////////////////////////////////////////////
  279.  
  280.             for (j=0; j<currentLength; j++) {
  281.                 negativeOneToI = -negativeOneToI;
  282.                 
  283.                 real_data[j] = vector_data[j] = 
  284.                     negativeOneToI + sin(3*PI*j/currentLength);
  285.             }
  286.  
  287.             /////////////////////////////////////////////////
  288.             // perform forward and inverse FFT and compare
  289.             // result with original signal
  290.             /////////////////////////////////////////////////
  291.  
  292.             result = FFTRealForward(vector_data, currentLength);
  293.             if (result != noErr) {
  294.                 printf("\nerror in forward real FFT");
  295.             }
  296.             
  297.             result = FFTRealInverse(vector_data, currentLength);
  298.             if (result != noErr) {
  299.                 printf("\nerror in forward real FFT");
  300.             }
  301.  
  302.             rmse_vector = RMSE_real(real_data, vector_data, currentLength);
  303.  
  304.             printf("\n real length %d FFT rmse: %g", currentLength, rmse_vector);
  305.         }
  306.         
  307.         if (real_data) DisposePtr((Ptr)real_data);
  308.         if (vector_data) DisposePtr((Ptr)vector_data);
  309.  
  310.     
  311.     }        
  312.  
  313.     
  314.  
  315. }
  316.  
  317. #define COMPLEX_CONVOLVE_MIN_POWER    2
  318. #define COMPLEX_CONVOLVE_MAX_POWER    13
  319. static void TestComplexConvolve()
  320. {
  321.     float     *data1;
  322.     float     *data2;
  323.     float    *literalConvolveData;
  324.     int        i, j;
  325.     double    negativeOneToI = -1;
  326.     double    rmse_vector;
  327.     long    currentLength;
  328.     OSErr    result;
  329.     
  330.     for (i=COMPLEX_CONVOLVE_MIN_POWER; i <=    COMPLEX_CONVOLVE_MAX_POWER; i++) {
  331.         currentLength = 1 << i;
  332.         
  333.         data1 =                 (float*)NewPtr(2*currentLength*sizeof(float));
  334.         data2 =                 (float*)NewPtr(2*currentLength*sizeof(float));
  335.         literalConvolveData =    (float*)NewPtr(2*currentLength*sizeof(float));
  336.         
  337.         if (data1 == nil || data2 == nil || literalConvolveData  == nil) {
  338.  
  339.             printf("\nerror allocating space for complex convolve test");
  340.  
  341.         } else {
  342.  
  343.             /////////////////////////////////////////////////
  344.             // initialize signals
  345.             /////////////////////////////////////////////////
  346.  
  347.             for (j=0; j<currentLength; j++) {
  348.                 data1[2*j]         =  cos(PI*j*j/currentLength);
  349.                 data1[2*j+1]     =  sin(PI*j*j/currentLength);
  350.                 
  351.                 data2[2*j]         =  2*sin(3.3*PI*j*j/currentLength);
  352.                 data2[2*j+1]     =  -cos(1.4*PI*j*j/currentLength);
  353.                 
  354.                     
  355.             }
  356.  
  357.  
  358.             /////////////////////////////////////////////////
  359.             // perform literal convolution, and compare
  360.             // with convolution performed through FFT
  361.             /////////////////////////////////////////////////
  362.  
  363.             ConvolveComplexLiteral(data1, data2, literalConvolveData, currentLength);
  364.  
  365.             result = ConvolveComplexAltivec(data1, data2, currentLength);
  366.             
  367.             if (result != noErr) {
  368.                 printf("\nerror in convolve complex");
  369.             }
  370.  
  371.             rmse_vector = RMSE_complex(data2, literalConvolveData, currentLength);
  372.  
  373.             printf("\n complex length %d convolve rmse: %g", currentLength, rmse_vector);
  374.         }
  375.         
  376.         if (data1) DisposePtr((Ptr)data1);
  377.         if (data2) DisposePtr((Ptr)data2);
  378.         if (literalConvolveData) DisposePtr((Ptr)literalConvolveData);
  379.  
  380.     
  381.     }
  382.     
  383.  
  384. }
  385.  
  386.  
  387.  
  388. #define REAL_CONVOLVE_MIN_POWER    0
  389. #define REAL_CONVOLVE_MAX_POWER    14
  390. static void TestRealConvolve()
  391. {
  392.     float     *data1;
  393.     float     *data2;
  394.     float    *literalConvolveData;
  395.     int        i, j;
  396.     double    negativeOneToI = -1;
  397.     double    rmse_vector;
  398.     long    currentLength;
  399.     OSErr    result;
  400.     
  401.     for (i=REAL_CONVOLVE_MIN_POWER; i <=    REAL_CONVOLVE_MAX_POWER; i++) {
  402.         currentLength = 1 << i;
  403.         
  404.         data1 =                 (float*)NewPtr(currentLength*sizeof(float));
  405.         data2 =                 (float*)NewPtr(currentLength*sizeof(float));
  406.         literalConvolveData =    (float*)NewPtr(currentLength*sizeof(float));
  407.         
  408.         if (data1 == nil || data2 == nil || literalConvolveData  == nil) {
  409.  
  410.             printf("\nerror allocating space for real convolve test");
  411.  
  412.         } else {
  413.  
  414.             /////////////////////////////////////////////////
  415.             // initialize signals
  416.             /////////////////////////////////////////////////
  417.  
  418.             for (j=0; j<currentLength; j++) {
  419.                 negativeOneToI = -negativeOneToI;
  420.                 
  421.                 data1[j] = 
  422.                     negativeOneToI + sin(3*PI*j/currentLength);
  423.                     
  424.                 data2[j] =     
  425.                     3 * negativeOneToI + sin(4.4*PI*j/currentLength);
  426.                 
  427.             }
  428.  
  429.             /////////////////////////////////////////////////
  430.             // perform literal convolution, and compare
  431.             // with convolution performed through FFT
  432.             /////////////////////////////////////////////////
  433.  
  434.             ConvolveRealLiteral(data1, data2, literalConvolveData, currentLength);
  435.  
  436.             result = ConvolveRealAltivec(data1, data2, currentLength);
  437.             
  438.             if (result != noErr) {
  439.                 printf("\nerror in convolve real");
  440.             }
  441.  
  442.             rmse_vector = RMSE_real(data2, literalConvolveData, currentLength);
  443.  
  444.             printf("\n real length %d convolve rmse: %g", currentLength, rmse_vector);
  445.         }
  446.         
  447.         if (data1) DisposePtr((Ptr)data1);
  448.         if (data2) DisposePtr((Ptr)data2);
  449.         if (literalConvolveData) DisposePtr((Ptr)literalConvolveData);
  450.  
  451.     
  452.     }
  453.     
  454.  
  455. }
  456.  
  457.  
  458. #define COMPLEX_FFT2D_MIN_POWER_X    1
  459. #define COMPLEX_FFT2D_MAX_POWER_X    8
  460. #define COMPLEX_FFT2D_MIN_POWER_Y    1
  461. #define COMPLEX_FFT2D_MAX_POWER_Y    8
  462.  
  463. static void TestFFT2DComplex()
  464. {
  465.     float     *data1;
  466.     float     *data2;
  467.     int        xsize, ysize, j;
  468.     double    negativeOneToI = -1;
  469.     double    rmse_vector;
  470.     long    currentLength;
  471.     OSErr    result;
  472.     
  473.     for (xsize=COMPLEX_FFT2D_MIN_POWER_X; xsize <=    COMPLEX_FFT2D_MAX_POWER_X; xsize++) {
  474.     
  475.         for (ysize = COMPLEX_FFT2D_MIN_POWER_Y; ysize <= COMPLEX_FFT2D_MAX_POWER_Y; ysize++) {
  476.         
  477.             currentLength = (1 << xsize) * (1 << ysize);
  478.             
  479.             data1 =                 (float*)NewPtr(2*currentLength*sizeof(float));
  480.             data2 =                 (float*)NewPtr(2*currentLength*sizeof(float));
  481.             
  482.             if (data1 == nil || data2 == nil ) {
  483.  
  484.                 printf("\nerror allocating space for FFT 2D test");
  485.  
  486.             } else {
  487.  
  488.             /////////////////////////////////////////////////
  489.             // initialize signal
  490.             /////////////////////////////////////////////////
  491.  
  492.                 for (j=0; j<currentLength; j++) {
  493.                     data1[2*j] = data2[2*j] =
  494.                         cos(PI*j*j/currentLength);
  495.                     data1[2*j+1] = data2[2*j+1] =
  496.                         sin(PI*j*j/currentLength);
  497.                         
  498.                 }
  499.  
  500.                 /////////////////////////////////////////////////
  501.                 // perform forward and inverse FFT and compare
  502.                 // result with original signal
  503.                 /////////////////////////////////////////////////
  504.  
  505.                 result = FFT2DComplex(data1, 1 << xsize, 1 << ysize, -1);            
  506.                 if (result != noErr) {
  507.                     printf("\nerror in fft 2d complex");
  508.                 }
  509.                 
  510.                 result = FFT2DComplex(data1, 1 << xsize, 1 << ysize, 1);
  511.                 if (result != noErr) {
  512.                     printf("\nerror in fft 2d complex");
  513.                 }
  514.                 
  515.                 rmse_vector = RMSE_complex(data1, data2, currentLength);
  516.                 
  517.                 printf("\n complex 2D %d X %d FFT rmse: %g", 1 << xsize, 1 << ysize, rmse_vector);
  518.                 
  519.             }
  520.             
  521.             if (data1) DisposePtr((Ptr)data1);
  522.             if (data2) DisposePtr((Ptr)data2);
  523.         
  524.         }
  525.         
  526.     }
  527.     
  528. }
  529.  
  530.  
  531. #define REAL_FFT2D_MIN_POWER_X    3
  532. #define REAL_FFT2D_MAX_POWER_X    9
  533. #define REAL_FFT2D_MIN_POWER_Y    2
  534. #define REAL_FFT2D_MAX_POWER_Y    9
  535. static void TestFFT2DReal()
  536. {
  537.     float     *data1;
  538.     float     *data2;
  539.     int        xsize, ysize, j;
  540.     double    negativeOneToI = -1;
  541.     double    rmse_vector;
  542.     long    currentLength;
  543.     OSErr    result;
  544.     
  545.     for (xsize=REAL_FFT2D_MIN_POWER_X; xsize <=    REAL_FFT2D_MAX_POWER_X; xsize++) {
  546.     
  547.         for (ysize = REAL_FFT2D_MIN_POWER_Y; ysize <= REAL_FFT2D_MAX_POWER_Y; ysize++) {
  548.         
  549.             currentLength = (1 << xsize) * (1 << ysize);
  550.             
  551.             data1 =                 (float*)NewPtr(2*currentLength*sizeof(float));
  552.             data2 =                 (float*)NewPtr(2*currentLength*sizeof(float));
  553.             
  554.             if (data1 == nil || data2 == nil ) {
  555.  
  556.                 printf("\nerror allocating space for FFT 2D test");
  557.  
  558.             } else {
  559.  
  560.                 /////////////////////////////////////////////////
  561.                 // initialize signal
  562.                 /////////////////////////////////////////////////
  563.  
  564.                 for (j=0; j<currentLength; j++) {
  565.                     negativeOneToI = -negativeOneToI;
  566.                     
  567.                     data1[j] = data2[j] = 
  568.                         negativeOneToI + sin(3*PI*j/currentLength);
  569.                 }
  570.  
  571.                 /////////////////////////////////////////////////
  572.                 // perform forward and inverse FFT and compare
  573.                 // result with original signal
  574.                 /////////////////////////////////////////////////
  575.  
  576.  
  577.                 result = FFT2DRealForward(data1, 1 << xsize, 1 << ysize);            
  578.                 if (result != noErr) {
  579.                     printf("\nerror in fft 2d real forward");
  580.                 }
  581.                 
  582.                 result = FFT2DRealInverse(data1, 1 << xsize, 1 << ysize);
  583.                 if (result != noErr) {
  584.                     printf("\nerror in fft 2d real inverse");
  585.                 }
  586.                 
  587.                 rmse_vector = RMSE_real(data1, data2, currentLength);
  588.                 
  589.                 printf("\n real 2D %d X %d FFT rmse: %g", 1 << xsize, 1 << ysize, rmse_vector);
  590.                 
  591.             }
  592.             
  593.             if (data1) DisposePtr((Ptr)data1);
  594.             if (data2) DisposePtr((Ptr)data2);
  595.         
  596.         }
  597.         
  598.     }
  599.     
  600. }
  601.  
  602. #define COMPLEX_CONVOLVE2D_MIN_POWER_X    2
  603. #define COMPLEX_CONVOLVE2D_MAX_POWER_X    5
  604. #define COMPLEX_CONVOLVE2D_MIN_POWER_Y    2
  605. #define COMPLEX_CONVOLVE2D_MAX_POWER_Y    5
  606.  
  607.  
  608. #define REAL_CONVOLVE2D_MIN_POWER_X    3
  609. #define REAL_CONVOLVE2D_MAX_POWER_X    8
  610. #define REAL_CONVOLVE2D_MIN_POWER_Y    2
  611. #define REAL_CONVOLVE2D_MAX_POWER_Y    8
  612.  
  613. static void TestConvolve2DReal()
  614. {
  615.     float     *data1;
  616.     float     *data2;
  617.     float    *literalConvolve;
  618.     int        xsize, ysize, j;
  619.     double    negativeOneToI = -1;
  620.     double    rmse_vector;
  621.     long    currentLength;
  622.     OSErr    result;
  623.     
  624.     for (xsize=REAL_CONVOLVE2D_MIN_POWER_X; xsize <=    REAL_CONVOLVE2D_MAX_POWER_X; xsize++) {
  625.     
  626.         for (ysize = REAL_CONVOLVE2D_MIN_POWER_Y; ysize <= REAL_CONVOLVE2D_MAX_POWER_Y; ysize++) {
  627.         
  628.             currentLength = (1 << xsize) * (1 << ysize);
  629.             
  630.             data1 =                 (float*)NewPtr(currentLength*sizeof(float));
  631.             data2 =                 (float*)NewPtr(currentLength*sizeof(float));
  632.             literalConvolve =        (float*)NewPtr(currentLength*sizeof(float));
  633.             
  634.             if (data1 == nil || data2 == nil || literalConvolve == nil) {
  635.  
  636.                 printf("\nerror allocating space for convolve 2D real test");
  637.  
  638.             } else {
  639.  
  640.                 /////////////////////////////////////////////////
  641.                 // initialize signals
  642.                 /////////////////////////////////////////////////
  643.  
  644.                 for (j=0; j<currentLength; j++) {
  645.                     negativeOneToI = -negativeOneToI;
  646.                     
  647.                     data1[j] = negativeOneToI + sin(0.05*PI*j/currentLength);
  648.                     data2[j] = -1.3  * cos(1.3*PI*j/currentLength);
  649.                 }
  650.  
  651.                 /////////////////////////////////////////////////
  652.                 // perform literal convolution, and compare
  653.                 // with convolution performed through FFT
  654.                 /////////////////////////////////////////////////
  655.  
  656.                 Convolve2DRealLiteral(data1, data2, literalConvolve, 1 << xsize, 1 << ysize);
  657.  
  658.                 result = ConvolveRealAltivec2D(data1, data2, 1 << xsize, 1 << ysize);
  659.                 if (result != noErr) {
  660.                     printf("\nerror in fft 2d real");
  661.                 }
  662.                 
  663.                 rmse_vector = RMSE_real(data2, literalConvolve, currentLength);
  664.                 
  665.                 printf("\n real 2D %d X %d convolve rmse: %g", 1 << xsize, 1 << ysize, rmse_vector);
  666.                 
  667.             }
  668.             
  669.             if (data1)                 DisposePtr((Ptr)data1);
  670.             if (data2)                 DisposePtr((Ptr)data2);
  671.             if (literalConvolve)     DisposePtr((Ptr)literalConvolve);
  672.         
  673.         }
  674.         
  675.     }
  676.     
  677. }
  678.  
  679.  
  680. #pragma mark -
  681.  
  682. int main() {
  683.     
  684.     ////////////////////////////////////////////////////////////////////////////////
  685.     
  686.     printf("\n\nTesting complex FFT:\n");
  687.     
  688.     TestComplexFFT(); 
  689.  
  690.     ////////////////////////////////////////////////////////////////////////////////
  691.  
  692.     printf("\n\nTesting real FFT:\n");
  693.  
  694.     TestRealFFT();
  695.  
  696.     ////////////////////////////////////////////////////////////////////////////////
  697.  
  698.     printf("\n\nTesting 2D complex FFT:\n");
  699.  
  700.     TestFFT2DComplex();
  701.  
  702.     ////////////////////////////////////////////////////////////////////////////////
  703.  
  704.     printf("\n\nTesting 2D real FFT:\n");
  705.  
  706.     TestFFT2DReal();
  707.  
  708.     ////////////////////////////////////////////////////////////////////////////////
  709.  
  710.     printf("\n\nTesting complex convolution:\n");
  711.  
  712.     TestComplexConvolve();
  713.  
  714.     ////////////////////////////////////////////////////////////////////////////////
  715.  
  716.     printf("\n\nTesting real convolution:\n");
  717.  
  718.     TestRealConvolve();
  719.  
  720.     ////////////////////////////////////////////////////////////////////////////////
  721.  
  722.     printf("\n\nTesting 2D real convolution:\n");
  723.  
  724.     TestConvolve2DReal();
  725. }
  726.